Two-dimensional Hubbard-Holstein bipolaron 
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We present a diagrammatic Monte Carlo study of the properties of the Hubbard-Holstein bipo- 
laron on a two-dimensional square lattice. With a small Coulomb repulsion, U, and with increasing 
electron-phonon interaction, and when reaching a value about two times smaller than the one cor- 
responding to the transition of light polaron to heavy polaron, the system suffers a sharp transition 
from a state formed by two weakly bound light polarons to a heavy, strongly bound on-site bipolaron. 
Aside from this rather conventional bipolaron a new bipolaron state is found for large U at interme- 
diate and large electron-phonon coupling, corresponding to two polarons bound on nearest-neighbor 
sites. We discuss both the properties of the different bipolaron states and the transition from one 
state to another. We present a phase diagram in parameter space defined by the electron-phonon 
coupling and U. Our numerical method does not use any artificial approximation and can be eas- 
ily modified to other bipolaron models with longer range electron-phonon and /or electron-electron 
interaction. 



I. INTRODUCTION 

The interaction between electrons and lattice degrees of 
freedom plays a crucial role in the properties of many ma- 
terials and results in a multitude of physical phenomena. 
Structural transitions like the cooperative Jahn- Teller 
distortion in perovskites, the Peierls dimerization in one- 
dimensional systems, or the pairing and the condensation 
of the charge carriers in superconductivity are some of the 
most spectacular effects which originate from electron- 
lattice interaction. Polaronic and bipolaronic effects are 
found in many materials like transition metal-oxidesi, su- 
perconducting materials^ and conjugated polymers*. In 
the last years there has been growing experimental ev- 
idence that even in the fashionable strongly correlated 
materials like manganates and cuprates, aside from the 
unscreened Coulomb repulsion the electron-lattice inter- 
action is extremely important^. 

The theoretical investigation of the interaction between 
charge carriers and lattice vibrations has a long history. 
The concept of a polaron, which describes an electron 
which carries with it a lattice deformation, was introduce 
first by Landau in 1933&. One of the most successful the- 
ories of the last century, the Bardeen, Cooper and Schri- 
effer (BCS) theory of superconductivity is based on the 
observation that phonons induce an effective attraction 
between electrons. 

The discovery of high T c superconductors renewed the 
interest in the study of the electron-phonon interaction, 
and the bipolaron problem in particular. In the cuprate 
materials the strength of the interaction between elec- 
trons and phonons is believed to be in the intermedi- 
ate regime. Because of this rather strong interaction the 
lattice ions change their equilibrium position when they 
are in the vicinity of charge carries, i.e the charge car- 
ries drive a phonon vacuum instability^*^. The classical 
Migdal-Eliashberg approach to the theory of supercon- 
ductivity, which neglects these effects and is valid only 



for small electron-phonon coupling, cannot be applied 
to cuprates. Therefore special theoretical attention was 
given to scenarios where electrons (or holes) form pairs in 
real space (bipolarons), which suffer Bose-Einstein con- 
densation, leading to superconductivity. In this respect 
Alexandrov and co-workers proposed the strong electron- 
phonon coupling theory as a starting point for explaining 
the physics of high T c superconductors. In their the- 
orji2*ii, as a consequence of the strong electron-phonon 
interaction, the holes pair in bipolarons. For low charge 
carriers density, the system can be regarded as a charged 
2e Bose gas which condense at T c , resulting in supercon- 
ductivity. 

However the study of these systems is complicated by 
the failure of both strong and weak coupling perturba- 
tion theory, even for simple model systems, at intermedi- 
ate coupling strength. Novel algorithms were developed 
to address the problem in this difficult region of param- 
eters space. The bipolaron problem, which is defined by 
two electrons on a lattice, has been intensively studied in 
the last years. For the one dimensional case the problem 
was addressed by Bonca et aL 13 with an exact diagonal- 
ization technique on a variationally determined Hilbcrt 
subspace. Other one-dimensional calculations were based 
on variational methods^ and density-matrix combined 
with Lancszos diagonalization techniquoi^. The two- 
dimensional case was investigated in the adiabatic ap- 
proximation 1 - 6 , and with variational methods^. 

In this paper we address the Hubbard-Holstein (HH) 
bipolaron which is one of the simplest and most pop- 
ular models which contains both electron-electron and 
the electron-phonon interaction. Its solution is important 
for understanding the competition between the phonon- 
induced electron-electron attraction and the electron- 
electron Coulomb repulsion 3 ^. As our calculation, based 
on a Quantum Monte Carlo algorithm, shows, in the in- 
termediate and strong electron-phonon coupling regions, 
phonons strongly renormalize the effective hopping in- 
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tegral of the electrons, strongly reduce the effective on- 
site Coulomb repulsion but do not significantly affect the 
nearest-neighbor exchange interaction. This gives rise to 
low-energy effective Hamiltonians with a large antifer- 
romagnetic interaction relative to the the effective hop- 
ping and the effective repulsion terms, which couldn't be 
derived starting from pure electronic models. We find 
that, depending on the value of the Coulomb repulsion, 
two electrons can form an on-site strongly bound state 
for small U or a weakly bound nearest-neighbor local- 
ized state for larger U. The former state appears when 
the effective on-site attraction due to phonons overcome 
the Coulomb repulsion and the later is a result of the 
exchange interaction which wins over the strongly renor- 
malized electron (polaron) kinetic energy. 

We developed a Diagrammatic Quantum Monte Carlo 
(DQMC) algorithm suitable for studying the two- 
dimensional HH bipolaron. To our knowledge, this is 
the first two-dimensional calculation which considers dy- 
namical phonons and does not entail any artificial trun- 
cation of the Hilbert space. Our algorithm computes the 
imaginary time two-particle Green's function from which 
we extract information about the bipolaron state at long 
imaginary time. The DQMC algorithms were introduced 
by Prokof 'ev et al*& and used to calculate the properties 
of Frohlich 19 and spin 20 polarons. With the same tech- 
nique the two-body problem was addressed by Burovski 
et alm^ for the exciton problem. Here we work in di- 
rect space (site) representation, the basis consisting of 
Wannier orbitals and phonons at each site. This is in 
contrast to the exciton problem where the electron-hole 
interaction is attractive allowing a momentum space cal- 
culation free of the sign problem. By working in real 
space we have managed to avoid the sign problem which 
would appear in the momentum representation when the 
Coulomb repulsion is introduced. The code can be eas- 
ily adapted to include longer range electron-phonon or 
electron-electron interactions and to study models more 
suitable to the cuprates, as for example the extended HH 
model 2 ^. The disadvantage is that in this basis the mo- 
mentum dependence of different quantities is difficult or 
sometimes impossible to compute. We consider a square 
lattice of 25 x 25 sites with periodic boundary condi- 
tions which is large enough for negligible finite size errors. 
There are no other truncations of the Hilbert space. 

II. MODEL HAMILTONIAN 

The Hubbard-Holstein Hamiltonian reads 

H = ~t ( C *W + H - c ) + U Y, "*T"i4. + 
(if) i 

+LO Q ^b\bi + g^n la {b\ + hi) . (1) 

i i,cr 

Here c i<7 (ci a ) is the creation (annihilation) operator of 



an electron with spin a at site i. b[, bi are phonon cre- 
ation and respectively annihilation operators at site i. 
The first term describes the nearest-neighbor hopping of 
the electrons, and the second the on-site Coulomb re- 
pulsion between two electrons. The lattice degrees of 
freedom are described by a set of independent oscillators 
at each site, with frequency ujq. The electrons couple 
through the density ni a = c\ a Ci a to the local lattice dis- 
placement Xi oc (pi + bi) with a strength g. This Hamil- 
tonian describes a tight-binding model together with an 
on-site Coulomb repulsion term and an on-site electron- 
phonon interaction term. The Holstein and the Hubbard 
models are limiting cases for [7 = and g = respec- 
tively. 

In this paper we address the electron pairing as a func- 
tion of both Coulomb repulsion and electron-phonon in- 
teraction by studying two electrons on a square two- 
dimensional lattice. 



III. PERTURBATION THEORY RESULTS 

The HH model cannot be solved analytically except 
for the two extreme cases of weak and strong electron- 
phonon interaction, when perturbation theory can be ap- 
plied. The most interesting physical situation is in be- 
tween these regimes. In order to understand what hap- 
pens in the intermediate region, it is necessary to present 
first the weak and strong coupling cases. 



A. Weak electron-phonon coupling 

For g = 0, the ground state will be formed by two 
electrons with zero momentum moving freely through the 
lattice. The total spin is zero because the triplet state 
could not have the two electrons in the same k = state. 

When the electron-phonon interaction is switched on, 
two things happen. The electrons get lightly dressed 
which results in an increase of their effective mass, and 
the electron-phonon interaction introduces a frequency 
dependent effective attraction between the electrons. Up 
to second order in g the effective attraction is propor- 
tional to the phonon propagator 

V e f f (u)= 9 >V( q ,u) = -^L . (2 ) 

This is a retarded interaction and attractive at small 
frequency (for to < loq). In the antiadiabatic limit 
(ujq — > oo) where the ions are considered light and able 
to follow instantaneously the motion of the electrons, the 
effective interaction (Eq. |5J) is instantaneous. 

In our model the Coulomb repulsion competes with the 
phonon-induced attraction, resulting in a total effective 
interaction 
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V eff (w) = U- 



2g 2 uj 



(3) 



In the antiadiabatic limit (loq — ► oo) the effective inter- 
action (Eq. 0) is instantaneous with 



eff 



U 



(4) 



and the situation can be described by a pure Hubbard 
model. 

The ability of the interaction to bind the electrons 
into a pair is essential. It is well known that in a two- 
dimensional lattice any attractive instantaneous interac- 
tion will cause two electrons to bind in a pair. We are 
not aware of any analytical proof that this is true for 
a retarded interaction (finite phonon frequency), but a 
simple numerical calculation on a variational space which 
allows states with maximum one phonon shows that the 
binding persists for Lu /t < 1. However, the binding en- 
ergy is always very small (smaller than 10~ 3 t) and de- 
creases roughly linearly with decreasing u>q. With in- 
creasing U, in antiadiabatic limit, the pairing persists 
as long as U < 2g 2 /too- For finite uio the binding en- 
ergy decreases rapidly with increasing U. It is plausible 
that the critical U defining the pair formation to be still 
2<7 2 /wo, but the binding energy becomes so small with 
increasing U that it is very difficult to resolve numeri- 
cally. However, we do believe that such extremely weak 
bound pairs have no physical importance. Besides, when 
the electron-phonon coupling is increased the situation 
gets rapidly complicated and the calculations restricted 
on the variational space with a maximum of one phonon 
are inappropriate. The renormalization of the electron- 
phonon interaction vertex becomes important. Migdal's 
theoremSi applied in classical superconductivity theory is 
not valid here because of the absence of the Fermi sea 31 . 



Cia=c ia e^ h \ b ^ . (8) 
The Hamiltonian written in the new basis is 

H = H t +H Q (9) 

with 



tf = Wo £^--E^>+( c/ -— )V 

* ' C/Jn ' C/Jn ' 



(10) 



H * = "* E (45^4^ + H -c) (11) 



and 



(12) 



The physical meaning of this canonical transformation 
is a shift of the ions equilibrium position at the sites 
where the electrons are present 



(xi> = (b\ + h) = (b\ +bi+J2 — «-) = te> + — (m) ■ 

(13) 

As can be seen from the second term of Eq.^JJ the lattice 
deformation energy gained due to the electron presence 
is 



E p = g 2 /uj 



(14) 



The dimensionless electron-phonon coupling constant 
may be defined as the ratio between this energy and the 
bare electron kinetic energy which is proportional to the 
hopping t and with the lattice dimensionality z. We de- 
fine it as2£ 



B. Strong electron-phonon coupling 

In the strong coupling limit, the HH model may be ad- 
dressed with a canonical transformation and by treating 
the hopping part of the Hamiltonian as a perturbation. 
The last three terms in Eq. are diagonal in the ro- 
tated basis obtained by applying the unitary operator— 

Q 

e where 



S 



9_ 

too 



^nur(bt-bi) 



Using formula 



A = e s Ae~ s = A + [S, A] + -[S, [S, A}} 

the transformed operators become 
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bi = bi 



(5) 



(6) 



(7) 



uiozt 2uiot 



(15) 



Since the electron hopping is accompanied by a change in 
ions equilibrium position (see the term XjXj in Ea. llljl . 
it is exponentially reduced for large g. 



t e f f = t(i\c\c j xjx j \j)=te k=te Z* (16) 



The effective on-site interaction between electrons is 

t 
to 



U, 



eff 



U-2 9 - 



U-2E, 



(17) 



(the same as in Eq. Q), and in the antiadiabatic limit 
(when Lu ,g — > cxd, g/u>Q — > and 2g 2 /uj^ is finite) 
Xi = I and the model can be mapped again in a pure 
Hubbard one. 
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FIG. 1: Schematic representation of second order processes 
in H t which determine the elements of matrix T. The small 
horizontal lines represent the lattice sites. If an electron is on 
a particular site then the equilibrium position of the ion at 
that site will be changed, and we indicate this by a downwards 
shift. The excitation energy of the intermediate state is shown 
above every process. If the initial and the final lattice config- 
urations are different, the corresponding matrix element will 
be exponentially reduced. In the large g limit the correspond- 

2 2/2 2 2/2 

ing matrix elements are: a) -j^~e~ 9 ' UJ °- b) -t^p~9 > u ° . 
c) - ,, 1 „ e 



E v 



For negative U e f f it is evident that the electrons form a 
bound state. However, based on second order perturba- 
tion theory in the hopping _ff t , it can be shown that even 
for positive U e ft a stable bipolaron state can exist. Let's 
consider the case of large U which results in U e ff > 0. 
The ground state of Hq is formed by the degenerate states 



k) = c iT c i+oil )' 



a f 



(18) 



The meaning of this notation is that the electron with 
spin J. is at a distance "a" from the one with spin \ re- 
siding at site "i" . "a" can take all the possible values 
except 0. In first order perturbation theory, the calcula- 
tion of the matrix (<Zj | H t \ bj) reduces to the calculation 
of Ea. 1 161 and indicates an exponentially reduced nearest- 
neighbor hopping. Second order perturbation theory sta- 
bilizes the bipolaron states. It is equivalent to diagonal- 
izing the operator 



T = H t 



1 



En — Hi 



H t 



(19) 



on the subspace spanned by all the degenerate states 
of Hq. The processes which can take place are shown 
schematically in Fig. ^ We can classify them in two 
classes. The first class includes processes as in Fig. ^-a, 
-b and -c, where the final lattice configuration is different 



from the initial one. This results in an exponential reduc- 
tion of the matrix elements, so that we can neglect them 
in first approximation. The second class includes pro- 
cesses like the ones in Fig. Q]-d and -e where the initial 
and the final lattice configuration is unchanged. They 
are not exponentially reduced. In Fig. ^-d an electron 
hops on a neighboring site without carrying the lattice 
deformation around it and afterwards comes back. The 
energy of the intermediate state is 2E p , because it con- 
tains a site with deformation and without electron, and a 
site with an electron and without deformation. The gain 
in energy is —t 2 /2E p . This process only contributes to 
the diagonal elements of T. In Fig.^-e one electron hops 
without carrying the lattice deformation on the neigh- 
boring site which is occupied by the other electron. The 
intermediate state contains a doubly occupied site which 
has a deformation corresponding to only one electron and 
a site with deformation and without electron, therefore 
the energy of this state is U e f f + 2E p — U. The final state 
can be identical with the initial state (Fig.[I}el), and the 
process contributes to the diagonal elements T$,s, or the 
final state can have the electrons interchanged (e.g.: from 
the initial |f, j> to the final j, |> - Fig. |I}e2), and this 
process contributes to the non-diagonal elements T$,-s- 
5 means nearest-neighbor here. The energy gain corre- 
sponding to each of the two processes shown in Fig. ^e 
is —t 2 /U. Neglecting the exponentially reduced terms in 
the calculation of T, what remains is the diagonal terms 
and the off-diagonal ones which connect the states | 8) 
and | —<5). 



t 

T a ,a = 8 X (-Sjr-) 

Ts - S = 6x (-^) + 2x (-Tj) 
T s ^s = 2 x 



a ^ 0,1 



5=1 



(20) 



Solving the secular equation, we find the condition for 
the bipolaron existence to be 



U < AE; 



p 



and the bipolaron binding energy 



4f 2 



A 6 = -— + — 



E n 



U 



(21) 



(22) 



Notice that even though the bipolaron exists up to a large 
value of U it is a weakly bound state when U e ff > 0. 

The physical interpretation of these results is simple. 
The energy given by Eq. |22 corresponds to a double de- 
generate singlet state formed by two electrons on nearest- 
neighbor sites. One state is a singlet along the X direc- 
tion and the other along the Y direction. In distinction 
to Hubbard model where the exchange energy can never 
win over the kinetic energy and therefore cannot bind two 
electrons, here the interaction with phonons results in a 
strong reduction of the electron bandwidth but not of the 
exchange energy because it implies virtual transitions of 
electrons on double occupied sites without carrying the 
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lattice deformation with them. Therefore now the ex- 
change energy can easily win and produce singlet bound 
states. However there is another effect which introduces 
an effective repulsion between two nearest-neighbor elec- 
trons and wins over the exchange energy when U > AE p . 
A virtual transition of an electron to an empty nearest- 
neighbor site without carrying the lattice deformation 
will lower its energy by —t 2 /2E p (Fig. [TJd). But if the 
nearest-neighbor site is occupied by the other electron 
this process is not possible resulting in an effective re- 
pulsion of t 2 /E p between two nearest- neighbor electrons. 
Therefore Eq. [221 reflects the competition between this 
effective repulsion and the exchange attraction equal to 
-U 2 /U. 

When the processes shown in Fig. 2] -a, -b and -c are 
taken into account the degeneracy of the two singlets 
is lifted, and two states are formed. This results in a 
ground state with s-wave (Ai g ) symmetry and another 
state with d-wave (-Big) symmetry. It should also be 
mentioned that if a positive next-nearest-neighbor hop- 
ping t' is introduced in the model, the d-wave symmetry 
state will be stabilized and it becomes the ground state 
when t' is large enough. 

Let's summarize the strong coupling regime physics, 
neglecting at the beginning the exponentially reduced 
terms. When U is small the ground state energy is 
U — 4E p and consists of two electrons located on the same 
site. The first excited state contains one more phonon 
and has an energy U — AE p + uiq (this is a N degenerate 
state because the phonon can be at any site). When U 
is increased and U — AE p + luq becomes larger than —2E p 
(which is the zero order energy of two electrons staying 
on different sites), the first excited state is a double de- 
generate nearest- neighbor singlet. When the hopping is 
switched on the low-energy physics can be described by 
the Hamiltonian 

H = -teff £ ( C U ff + ff.c) + J$> Sj -^) 

+ V^mrij + U ef fY^ n it n il + H ' ( 23 ) 

(id) 1 

2 / 2 

with the hopping t e ft — t e J ' 0, the exchange J = 

At 2 t 2 

, the nearest- neighbor repulsion V — and the on- 
site interaction U e ff = U — 2E P . H' describes the pro- 
cesses shown in Fig.^ -0 ! -b and -c. Their magnitude is 

either or jjj^^ (see the caption of Fig. [I}, which 
is much smaller than t e f / . In the literatur e) 13 ! 16 the bipo- 



laron with the electrons located on the same site is called 
50, and the one with the electrons located on nearest- 
neighbor sites 51. The Hamiltonian ill'MI) describes the 
transition from the SO to the SI bipolaron in strong cou- 
pling regime. For small (i.e. negative) U e ff the order 
of the lowest energy states is: s-wave SO, s-wave SI, d- 
wave SI. When U e ff increases the SO state starts mixing 
with s-wave 51 state. The mixing between the 50 and 
the s-wave 51 states is of order of t e ff, and the splitting 
between the s-wave and the d-wave 51 states is given by 
H', thus being much smaller. The order of low-energy 
states becomes: linear combination of s-wave 50 and 51, 
d-wave 51, linear combination of s-wave 50 and 51. For 
larger U only two bound states exists: s-wave 51 and 
d-wave 51. In conclusion, the ground state evolves ana- 
lytically (crossover) from 50 bipolaron to s-wave 51 bipo- 
laron with increasing U. The situation is different for the 
first excited state. Here at a critical value of U, a nonan- 
alytical transition takes place, and the first excited state 
changes from s-wave symmetry to d-wave symmetry. 



IV. ALGORITHM 

A. General approach 

Our algorithm calculates different imaginary time 
Green's functions and relies upon the ability to project 
out the ground state properties by extrapolating to long 
complex times. Let's consider the equation 

^\e- TH \^)=Y,Mv)\ 2 e- T£ » (24) 

V 

where l^) is a whatever state and {|^)} form the complete 
set of the eigenstates with energies £ v . We see that at 
large r Eq. [31 converges to 

|(i/o |V) I 2 e- T ^° (25) 

where \vq) is the ground state of the system. Suppose 
the ground state is separated from the first excited state 
by a gap A. We can obtain the ground state energy and 
the overlap of the ground state with with an accuracy 
better than 1% (for example) calculating Eq.[2]at a time 
t w 5/A. 

Because the total momentum K is a quantity which is 
conserved in our problem we can obtain the lowest energy 
in the K channel by calculating 



pn ( K , T ) = £ (( K ~ k ~ qi ~ ■■■ ~ 1n)i,kf,qi, -,q n \e tH \(K - k - qi - ... - q n ) l ,k- l ;q 1 , ...,q n ) 

k,qx,...q n 

— > ]T \((K-k- qi - ...-q n ) l ,k r ,q 1 ,...,q n \v QK )\*e- TE W 

k.qi...q n 
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(26) 



r 



at large r. Here \ki±, k-zy, qi, q%, q n ) is a state with two 
electrons, one with momentum k\ and spin down and 
the other with momentum k% and spin up, and with n 
phonons with momentum q%, q2,--- and q n respectively. 
\vok) is the ground state (the state with the lowest en- 
ergy) in the K channel. The calculation of P n (K,r) 
yields both the bipolaron energy and the n-phonon con- 
figuration probability in the bipolaron state. 

For reasons related with the sign problem (to be dis- 
cussed later), we calculate P n (K,r) in real-space repre- 
sentation 



P "^' T ) = ^ E e iKx (i\e- TH T x \i) (27) 



i,X,ll ,1*2, ...l n 



where 



\i) = \i u (i + a)f,i + h,i + l 2 ,...,i + l„) (28) 



is a state with a spin down electron at site i, a spin up 
electron at site i + a and phonons at sites i + h, i + fa, 
. . . and i + l n , and 

T x \i) = 

\i + xi, (i + x + <x)t;« + x + h, i + x + fa, i + x + l n ) 

(29) 

is the state \i) translated with the vector x. 

Another quantity which is conserved is the total spin. 
Therefore for the singlet we calculate 

P?(K,t) = ± J2 e iKx {i s \e^ H T x \i s ) (30) 



with 



\i s ) = + a) s ;i + li,i + fa,...,i + l n ) (31) 



where (i,i + a) s is the singlet state with electrons at sites 
i and i + a. Similar equations can be written for the 
triplet channel. 

In order to calculate the above quantities, we developed 
a DQMC code which stochastically generates terms of the 
form 



G ij (r) = (i\e^ H \j(i,x)) 



(32) 



where \i) is a general state as in Eq. EHl with two 
electrons at arbitrary distance from each other and 
with an arbitrary number of phonons. The state 
\j(i,x)) can be obtain either by applying a transla- 
tion operation with an arbitrary vector x on \i) (i.e. 



\j) = T x \%i, (i + a) -[•; phonons)) or by applying a permu- 
tation (interchanging the electrons position) and a trans- 
lation on i) (i.e. \j) = T x \(i + a)j_, (i)^; phonons)). 

The value of an observable A in a particular K and S 
channel is 



e lKx ( m *> g s ' (m)w(m)a(m) 
Em w ( m ) 



(33) 



In Eq. 1331 we sum over all m generated configurations 
with the weight w(m) . M is the total number of measure- 
ments, x(m) is the translation vector which correspond 
to the configuration m, a(m) is the estimator of A and 
g s (m) is the factor which separates the triplet from the 
singlet. For singlet g (m) is 1 when electrons are on the 
same site and 1/2 otherwise. For triplet g s (m) is zero 
when the electrons are on the same site, and otherwise, 
1/2 when |j) is a translation of \i) and —1/2 when |j) is a 
translation of a state obtained from \i) by interchanging 
the electrons position. 



B. Implementation 

Let's start from the Hamiltonian QJ, and consider 



(34) 



as the noninteracting part of the Hamiltonian. Hq is 
diagonal in the real space representation. The evolution 
operator can be written as 



e-™ = e- THo S{T) 



(35) 



with 



( —\\ n r T r T 
S{ - T ) = Y.—T - dn...dT n T[H 1 { n )...H 1 {T n )} 
n=0 n\ Jq Jo 



(36) 



where Hi = H — Hq in the interaction picture is 
fli(r) ^e THo H 1 e- rHo . 

Eq. |22 becomes: 



(37) 
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Gij(r) = (*\e- TH \j) = e-^ £ [ —L / / ... / dn...dr n T[(i\H i (T 1 )H 1 (T 2 )...H i (r n )\j)] 

r, n - JO JO JO 



(38) 



r 
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FIG. 2: a) A typical diagram which represents a term in 
Eq. 13811 . The solid line (dashed line) represents a spin up 
(down) electron. The wavy line is a phonon propagator. The 
end at time r is a translation with the vector x of the end at 
time 0. b) The final state is a translation of the state obtained 
from initial state with the electrons position interchanged. In 
the code we consider only configurations where the final state 
is a translation of the initial state as in a) or is a transla- 
tion of the state obtained from initial state with the electrons 
position interchanged as in b). 
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this diagram is 
T3 ) The weight of 



a diagram is determined by the following rules: (i) every 
electron hopping corresponds to a term t dr, (ii) every 
electron-phonon vertex corresponds to a term g dr, (iii) 
every phonon propagator of length At corresponds to a term 
e -woAr^ anc j j-j v ~j evel y interval where the electrons are 
on the same site during a time At corresponds to a term 
-UAt 



the diagram ends must be either identical (Fig. |2J-a) or 
with the electrons position interchanged (Fig.|3-b). The 
rules for determining the diagram weight are given in the 
caption of Fig. [21 

We generate all the possible diagrams with an arbitrary 
number of phonons, and with the difference between 
and a T max chosen large enough to project the ground 
state. Our code is a continuous time code (i.e. it does 
not require artificial discretization of the imaginary time 
axis), and the diagrams are generated is a manner similar 
to that described in 18 i 19 i 21 . 

Estimators for energy, effective mass, phonon distri- 
bution and correlation function of the electrons position 
can be easily found. The measurements are taken only at 
large time where the ground state is projected out. The 
bipolaron energy estimator isi£ 



N„ 



-^V hop) 

(39) 

where r m is the time (length) of the m diagram, i v h 
counts the phonons propagators with Tj h length, jjj 
counts the intervals with double occupied sites with Tj u 
length, N ver tex is the number of electron-phonon vertices 
and Nhop is the number of electron hopping jumps. The 
estimator for the inverse of the bipolaron effective mass 
is 



2m„ 



(rn) 



x(mY 



(40) 



where x{m) is the translation vector between the time 
end and the time r m end. m e is the free electron effective 
mass. The probability to have n phonons is calculated 
with the estimator 



z n (m) 



(41) 



where n m is the number of phonons at the ends of the di- 
agram. The electrons relative positions correlation func- 
tion defined as 



CO) = Jj^{ n t n i+r) 



(42) 



has the estimator 



The calculation of Gij (r) is reduced to a series of in- 
tegrals, with an ever increasing number of integration 
variables. It is easy to show that every term in Eq. I|38|) 
can be represented by a diagram and a set of simple rules 
can be derived to determine the diagrams weight. Typ- 
ical examples of such diagrams are presented in Fig. |2 
Aside from a translation, the electronic configuration at 



C(r;m) = 6 Ti , 



(43) 



where R m is the relative distance between electrons at 
the ends of the measured diagram. 

As an illustration of our method, in Fig. ^ we show 
Iti(P s (0,t) * e^ T ) versus r where 



p s (o,T)*e^ = ]Tp;(o,T) 



(44) 
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FIG. 4: ln(P B (0,T) * e^ T ) versus r. a) For the given pa- 
rameters a large bipolaron binding energy results, b) For the 
given parameters a small bipolaron binding energy results. 
The linear asymptotic behavior starts at smaller time in case 
a). Notice that the time scale is different for the two cases 
presented. 

with P™ defined in Eq. for the K = channel, fi 
is an arbitrary parameter which is chosen close to the 
bipolaron energy to avoid the exponentially small weight 
of the large time diagrams. It can be seen that at long 
imaginary time ln(P s (0, r)) becomes linear in r, the slope 
being proportional to the ground state energy. An impor- 
tant remark should be made about the strong drop seen 
in P s (0,t) at short time. This is due to the fact that 
we generate only connected diagrams (diagrams where 
the phonon propagators are always glued to the electron 
propagator). The disconnected diagrams have an expo- 
nentially small contribution at large time, therefore they 
can be safely neglected, but at small time their omission 
will result in a strong potential drop which will not al- 
low an efficient sampling for both long and short time 
diagrams. We have eliminated this problem using a fic- 
titious potential renormalizatior. 18 . 



C. Discussions 

In order to avoid the sign problem, our code calcu- 
lates the Green's functions in the real-space representa- 
tion (Wannier basis) , even though this representation has 
some disadvantages. A similar expression to Eq. 1381 can 
be written in the momentum space (Bloch basis), and 
similar rules for determining the diagrams weight can be 
found. This approach was considered inSi for the exciton 
model calculation. The problem with our model is that, 
unlike in the exciton case where the conduction-electron 
valence-hole interaction is attractive, we have a repulsive 
interaction. This will make all the diagrams with an odd 
number of electron-electron interaction vertices negative, 
which implies a very severe sign problem. In the real- 



space representation the sign problem is avoided, all the 
diagrams being positive definite. However, this represen- 
tation introduces other problems, for example it makes 
the study of the bipolaron at large momentum difficult 
and inaccurate. In a more general sense, problems appear 
in all the irreducible channels aside from the one which 
contains the system ground state (i.e. K = 0, singlet). 
The difficulty is two-fold. First there is the sign problem. 
Aside from the singlet and K = channel where all the 
terms in Eq. 1331 are positive definite, at if ^ or/and in 
the triplet channel the factors e iKx(m) an( ^ g s (m) can 
take negative values. Second, in the real-space represen- 
tation we generate all the possible configurations with all 
the possible symmetries. To project out the lowest en- 
ergy state of the channel "7" we have to calculate P 7 (r) 
up to a time proportional to 1/A 7 , where A 7 is the 7 
channel gap. Therefore we have to simulate larger imag- 
inary times for a channel characterized by a smaller gap. 
But because all the symmetry channels configurations are 
generated in the same run, the statistics for the channel 
7 is proportional to e^^i ~ ^o) T (here Eq is the ground 
state energy). Thus if the imaginary time is increased, 
the statistics for channels other than the one which con- 
tains the ground state will be exponentially reduced. 

Another problem, specific to all ground state projec- 
tive algorithms, occurs if there are more than one bound 
state in the same symmetry channel, quasi-degenerate in 
energy. In this case, at large imaginary time we project 
out all these states. The results obtained in this case are 
going to reflect the average of the properties of all the 
projected states. We encounter this problem at large 
electron-phonon coupling, where the difference of the 
s-wave and c?-wave bipolaron energies is exponentially 
small, as the strong coupling theory predicts. However 
in the intermediate coupling regime we managed to sep- 
arate these states and we always found a s-wave ground 
state. 

Before we discuss the necessary modifications for 
adapting our algorithm to other bipolaron models, we 
mention that, even for the present HH bipolaron, the 
code can be improved. The momentum K and the spin 
S are not all the quantum labels which can be used to 
distinguish between the different symmetry channels. We 
also have the point group symmetries which break the 
Hilbert space in different irreducible representations. We 
have already discussed about s-wave {Ai g ) and rf-wave 
{Big) bipolaron states. In principle we can look for all 
the symmetries given by the representations of the D^h 
point group. These symmetries can be separated in a 
similar way to what we did when separating the singlet 
and the triplet channels. We have to generate diagrams 
where the time r is obtained after a translation and a 
point group operation applied to the time 0. In other 
words the electronic configurations of the diagrams ends 
should be connected by a space group operation^. For 
every operation there is a certain factor t D (m) which sep- 
arates the different representations, and the value of an 



observable is calculated analogues to Eq. 

A(K) = ^J2e iKx ^g s {m)t D (m)a(m) = 

rn 

For example, at K — 0, t D is always 1 for Ai g represen- 
tation. In general t D should be proportional to the char- 
acteristic of the representation. The sign problem can 
intervene for other representations than A\ g . The above 
approach applied to our present model will improve the 
accuracy of the results which describe the properties of 
the s-wave SI bipolaron. However the study of the d- 
wave symmetry bipolaron will still be difficult and not 
very accurate, because of the smallness of the binding 
energy of this state. We believe that no new physics will 
appear to justify the large coding effort necessary for im- 
plementing the above approach. Nevertheless the separa- 
tion of the different point group representations could be 
essential for other models which include longer electron- 
phonon interaction and where strongly bound bipolarons 
with electrons residing on different sites exisA 

The code can be easily modified to include longer range 
electron-electron and electron-phonon interaction. For 
an electron-phonon interaction term of the form 

Y^diM^ + bj) (46) 

and for an electron-electron interaction of the form 

X! 1 '</"'"• ( 47 ) 

there will be no sign problem. The diagrams are similar 
to the ones shown in Fig.|3aside from the possibility that 
now a phonon propagator at site "j" can be created or 
destroyed by an electron at site "i". This process will 
have a probability equal to <?y dr. If the two electrons 
are on sites "£" and "j" for an interval of time r a term 
e ~VijT j n ^ ne (jj a g ram W eight will correspond. 

V. RESULTS 

In all the subsequent calculations the phonon energy 
and the electron hopping term are chosen to be one (loq = 
1, t = 1). In real materials loq is smaller than t (by almost 
one order of magnitude), but calculations with such small 
wo are very time consuming and we find that they do not 
provide any new qualitative results. Our calculation is 
done on a 25 x 25 square lattice with periodic boundary 
conditions. 



A. Phonon Induced Attraction. U — Case 

With [7 = the effective interaction is always attrac- 
tive. In Fig. [S] we show the evolution of the system from 
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FIG. 5: U = 0, uio = 1) t = 1. a) The bipolaron energy versus 
electron-phonon coupling (circles). The dashed line (squares) 
is 2x free polaron energy, b) The bipolaron average number 
of phonons (circles) and 2 x free polaron number of phonons 
(squares), c) The probability to have the electrons on the 
same site, C(0), in the bipolaron state. The dotted line in -a 
( -b) represents the energy (number of phonons) of two free 
polarons versus the effective coupling, a e ff = 2 a. 




FIG. 6: A comparison between HH bipolaron (solid line) and 
attractive Hubbard model (dashed line). On the horizontal 
axis is U e ff/t e ff, where U e ff = — 4at and t e ff — t e ff(a) is 
the polaron effective hopping, a) The binding energy in terms 
of t e ff. b) Probability to have two electrons on the same site. 

a very weakly bond bipolaron (almost two free polarons) 
to a strongly bound bipolaron with increasing a. The 
transition is sharp. It can be seen from Fig. [5J-C, where 
the electrons relative position correlation function C{r) 
(Eq. 142(1 is shown at r = 0, that in a very short in- 
terval around a c = 1 the system ground state changes 
from almost two free polarons to a state where the elec- 
trons are practically on the same site. In comparison to 
the polaron case, the critical electron-phonon coupling 
where the transition takes place is approximately two 
times smaller. This can be understood by noticing that 
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(see Ea. llOf) the deformation energy at a particular site is 
proportional to the square of the number of electrons on 
that site. Therefore for the on-site bipolaron the effective 
a is two times larger than the corresponding polaron one. 
Thus, the bipolaron energy at a particular coupling a (in 
the strong coupling regime) is equal to two times the free 
polaron energy corresponding to a double a. The same 
is true for the average number of phonons in the bipola- 
ronic cloud (which is also proportional to the square of 
the number of electrons). This is shown with dotted line 
in Fig. -a and -b. From the same figures it can also 
be observed that the bipolaron transition is very sharp 
compared to the large to small polaron transition—. 

In Sec. IIIII we have shown that in the antiadiabatic 
limit (i.e when lu — ► oo) the effective attraction in- 
duced by phonons becomes instantaneous and as a con- 
sequence the HH model is equivalent to a an attractive 
Hubbard model. The attractive Hubbard Hamiltonian 
was considered as a realistic model to explain the proper- 
ties of systems like amorphous semiconductors^SSiSi or 
high T c superconductors^ and was under investigation in 
the pastS^. However, when the phonon frequency is finite, 
the interaction becomes retarded and the HH physics will 
differ from the attractive Hubbard one. Therefore, we 
think that a comparison of the HH model and the at- 
tractive Hubbard model is necessary. Aside from induc- 
ing a retarded effective attraction, the other main effect 
of the electron-phonon interaction is to dress the elec- 
trons, increasing their effective mass or equivalently re- 
ducing their effective hopping. Consequently, in Fig. 
we compare the HH Hamiltonian (with U — 0) with 
the corresponding attractive Hubbard model, defined by 
U e ff = —Aat = -2E p (see Eq. 0]and Eq. [TJ3 and po- 
laron effective hopping t e ff — t e ff(a). The effective 
polaron hopping as a function of electron-phonon cou- 
pling was calculated numerically with a DQMC code. In- 
creasing U e f f /t e f f the system evolves in both cases from 
two free electrons state to a state where the electrons 
are mainly on the same site (50 bipolaron). At small 
U e f f /tef f ' namely when a is smaller than the transition 
critical a Cl the effective attraction induced by phonons is 
weaker than the corresponding Hubbard attraction. The 
transition to the SO bipolaron is very sharp for the HH 
model, unlike the attractive Hubbard model where it is 
rather smooth. We found (not shown) that the smaller 
ujo, the sharper the transition is. This is also in agree- 
ment with the adiabatic limit calculation (u>o/t — > 0) 
where the transition is of first order—. Not only is the 
transition different but also the properties of the bipo- 
laron at large coupling are different in the two cases. The 
HH SO bipolaron has a large effective mass proportional 

to t^jj , thus with a factor of e^ a l u>0 larger than the free 
electron massif. On the other hand, the effective mass of 
the attractive Hubbard bipolaron increases linearly with 
|[/ e //| = 4at in the large U regime, as it can easily be 
shown analytically. 

From both Fig. Eland Fig. El we can conclude that for 
the HH model there is a very narrow transition region 



where the system evolves from two almost free light po- 
larons to a very heavy SO bipolaron. Before the tran- 
sition, the effect of the electron-phonon interaction is 
small, especially when the phonon frequency is small. In- 
creasing the phonon frequency results in increasing effec- 
tive attraction. The physics after the transition is well 
described by the strong coupling theory. Now the en- 
ergy and the number of phonons is proportional to a 
(E = —8at and N p h = 8at/uj ), as it can be seen in 
Fig. 0-a and -b, and the effective mass is proportional to 
e 4a/u;o u n lik e the weak coupling regime, here a smaller 
wq results in a heavier bipolaron state. 



B. SO Bipolaron to SI Bipolaron Transition. U ^ 

Case 

The weak coupling regime is characterized by the for- 
mation of a weakly bound state. The binding energy is 
extremely small even for U = 0, as can be seen from 
Fig. As discussed in Sec. IIII Al the bipolaron bind- 
ing energy decreases rapidly with increasing U, being 
well bellow the resolution limit of our algorithm (10 -3 t). 
Therefore we are not able to determine the critical U 
where the binding energy reaches zero. However, we do 
not consider this to be a relevant problem, a bipolaron 
with such a small binding energy being physically iden- 
tical with a state of two free polarons. 

In the strong coupling regime, as discussed in Sec. IIII Bl 
with increasing U the system evolves from a strongly 
bound SO to a weakly bound s-wave 51 bipolaron. The 
S0-S1 transition takes place around U = 2E p . At the 
critical value U = AE p the bipolaron state ceased to ex- 
ist, and the system becomes two polarons moving freely 
on the lattice. The binding energy of S'O bipolaron de- 
creases linearly with U. The SI bipolaron results from 
the exchange process and its binding energy is propor- 
tional to 1/17 fEa. 122)1. 

The energy of bipolaron in the intermediate coupling 
regime is presented in Fig. [7] When U = the bipolaron 
is a SO state. With increasing U the binding energy at 
first decreases linearly with U and afterwards its behav- 
ior changes, decreasing much slower. The bipolaron state 
disappears well before U reaches 4E P . The proportional- 
ity to U of the binding energy is a characteristic of the 
strongly bound S'O bipolaron. From Fig- [SI it can be seen 
that this is correlated with the probability to have the 
electrons on the same site, C(0) being close to one. At the 
value of U where the binding energy behavior changes, 
there is a transition to a weakly bound state (with the 
binding energy smaller than the one given by the strong 
coupling theory, Eq. where the probability to have 
the electrons on neighboring sites is enhanced and where 
simultaneously C(0) drops to small value. For large cou- 
plings, like in Fig.lSJa, this state is a small SI bipolaron, 
with the electrons residing essentially only on the nearest- 
neighbor sites. For smaller couplings, as in Fig. [SJ-d, this 
state is a "large" SI bipolaron, the wave function being 
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FIG. 7: a) The bipolaron energy versus U for different values 
of the electron-phonon coupling a. b), c) and d) The binding 
energy of the bipolaron in the transition region, for a = 1.62, 
a = 2.42 and respectively a = 3.125, versus U. The dashed 
line is the SI strong coupling theory predicted binding energy. 
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FIG. 8: The electrons position correlation function C(0) , 
4 x C(l), and 4 x C(2) (see Eq. versus U for different 

values of the electron-phonon coupling. 



spread over many sites, but still with an enhanced prob- 
ability that the electrons are nearest-neighbors. This can 
be seen from Fig.|3 where the correlation function, C(r) 
(Ea. I42|) . as a function of the electrons relative distance, 
r, is shown^. 

We know that, in the strong coupling regime for large 
U, a d-wave SI bipolaron, quasi-degenerate in energy 
with the s-wave SI ground state exists. When a is large 
(e.g. a = 3.125) our code projects out both SI states 
at large imaginary time. The results we obtain in this 
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FIG. 9: The correlation function C(r) in the intermediate 
electron-phonon coupling regime. The relative distance be- 
tween electrons is given in circular coordinates. 



case represent the average of the corresponding two SI 
bipolarons properties. At smaller couplings the energy 
difference between the two SI bipolarons is larger and 
it is easier to project out the ground state and thus to 
separate the two states. We show this in Fig. ^JJ where 
a comparison of the correlation function C{r) measured 
at two imaginary times, r = 35 and r = 80 is made. At 
time r = 35 we see a smaller probability for the electrons 
to stay along the diagonal directions. This is evidence 
that our measurements capture both the d-wave and the 
s-wave states. When the measurements are taken at a 
larger time, the value of the correlation function at sites 
which correspond to the diagonal directions increases. 
For the presented case, C(r) does not change sensible if 
the measurement time is increased above r = 50, thus we 
can conclude that the asymptotic regime is reached above 
this time. The fact that at r — 35 we see a decrease of 
C(5), i.e. a decrease of the correlation function at large 
distance along the diagonal directions, shows that the d- 
wave bipolaron in the intermediate coupling region is a 
large state spread over many sites, like the s-wave ground 
state. 

In Fig. ^2 we address the .S'O-.S'l transition by looking 
at the number of phonons in the bipolaron cloud and at 
bipolaron effective mass. Before the .SO-Sl transition, 
the number of phonons decreases slowly with U, and it is 
well approximated by the strong coupling perturbation 
theory. After the transition the number of phonons is 
roughly equal to the number of phonons corresponding to 
two free polarons. The S'O-ST transition is sharp, and the 
bipolaron changes from a very heavy state (5*0) to a light 
one (<S1), as can be seen in Fig.^^b where the inverse of 
the bipolaron effective mass is plotted as a function of U . 
The SQ-S1 transition is sharp even in the intermediate 
coupling region. We also have found (not shown) that 
for smaller uq the transition is sharper. However, we 
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FIG. 10: The correlation function C(r) in the intermediate 
electron-phonon coupling regime for a value of U which cor- 
responds to a large SI bipolaron. The sites at r — 2, r — 5 
and r — 9 are along the diagonal directions. The statistics 
are collected at two different imaginary times. For the smaller 
t — 35 time the results represent the average of the s-wave 
symmetry ground state and the d-wave symmetry first ex- 
cited state (notice the small occupation probability of the sites 
along the diagonal directions). At the larger time, r = 80, 
only the s-wave ground state remains. 
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FIG. 12: Phase diagram. The solid line is the bipolaron -two 
free polarons boundary. The dashed line separates the SO and 
SI bipolarons. The doted lines are the strong coupling theory 
results. 



smaller binding energy. The spatial extent of this state 
is also large. When a is increased both s-wave SI and 
d-wave ST bipolaron wave functions become more local- 
ized evolving to the states predicted by the perturbation 
theory. The energy difference between the s-wave and 
the d-wave states becomes exponentially small at large 
electron-phonon coupling. 

The results presented up to here are calculated for 
K = and in the singlet channel. For the strongly bound 
50 bipolaron, in agreement with the exponentially small 
effective hopping predicted by the perturbation theory, 
we found a flat dispersion resulting in a narrow band, 
of the order of the calculation error bars. We are not 
able to compute the momentum dependent properties 
of the weakly bound bipolarons for reasons described in 
Sec. IIV CI In the triplet channel we found no bound 
states for any value of the parameters. 



FIG. 11: a) Average number of phonons in the bipolaronic 
cloud for different values of electron-phonon coupling versus 
U. b) Inverse of the bipolaron effective mass versus U. m e 
is the electron mass. The horizontal lines correspond to the 
inverse of two free polarons effective mass. 

want to specify that when we talk about transition we 
mean a continuous change of system properties and not 
a non-analytical jump. 

To conclude, in the intermediate coupling regime, with 
increasing U, the system evolves continuously from a 
heavy, strongly bound 5*0 bipolaron to a light, weakly 
bound state spread over many sites. This state, which 
we call large SI bipolaron, has s-wave symmetry and an 
enhanced probability that the electrons occupy nearest- 
neighbor sites. In the same region of parameters an- 
other stable state with cZ-wave symmetry exists, with a 



C. Phase Diagram 

The phase diagram is shown in Fig.^| We want to re- 
mind the reader that with our technique, the calculation 
becomes difficult when the binding energy is small. The 
smaller is the binding energy, the larger time computa- 
tions are needed. The most difficult computations are 
at both large electron-phonon coupling and large U. The 
large imaginary time simulations are difficult because the 
number of phonons and the effective mass is always large 
in this case (the ground state consists of two weakly in- 
teracting small polarons, and a small polaron has itself 
an exponentially large mass and contains a large num- 
ber of phonons). Therefore the largest errors we get are 
in the determination of the bipolaron-two free polarons 
boundary at large values of a. In the strong coupling 
theory this boundary is determined by the critical value 
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U = 8at = 4Ep. In the intermediate coupling regime 
we find that the bipolaron state disappears much before 
that value. The value U = 8at = AE p should be taken 
as an upper limit for the existence of the bipolaron state, 
reached asymptotically when a is increased. 

The SI region contains a weakly bound state where 
the probability to have the electrons on nearest-neighbor 
sites is larger than of having them on the same site. De- 
pending on the value of electron-phonon coupling the SI 
bipolaron can be a large state with the wave function 
spread over may sites, or a small state where the elec- 
trons are residing on nearest-neighbor sites. The large 
SI bipolaron breaks into two large or intermediate (tran- 
sition) polarons and the small SI bipolaron evolves into 
two small polarons with increasing the Coulomb repul- 
sion U. The SI bipolaron is a state which forms only 
at intermediate and large electron-phonon coupling, thus 
only where the polaron kinetic energy is strongly reduced. 
In this region, the exchange attraction which is weakly 
renormalized by the phonons, can overcome the effective 
kinetic energy resulting in the formation of the SI bound 
state. 

At small a, the binding energy is extremely small, well 
beyond the resolution limit of our algorithm, and de- 
creases rapidly with increasing U . The maximum crit- 
ical U where the bipolaron breaks can be theoretically 
as large as U = 4at = 2E p , but for U larger than the 
boundary shown in Fig. ll2l thc binding energy is so small 
that can be safely approximated by zero. 

VI. CONCLUSIONS 

In this paper we studied the two-dimensional HH bipo- 
laron using a Diagrammatic Quantum Monte Carlo al- 
gorithm which computes the zero temperature Matsub- 
ara Green's functions. The bipolaron properties are 
extracted from the Green's functions behavior at large 
imaginary time where the ground state is projected out. 
Unlike the other DQMC simulations used for study- 
ing different polaron and exciton models, in order to 
avoid the sign problem, our algorithm produces and sums 
real space (Wannier orbitals basis) diagrams. The code 
can be relative easily modified for other bipolaron mod- 
els with longer range electron-phonon and/or electron- 
electron interaction. The dimensionality and the lattice 
symmetry can also be modified. 

We calculated the phase diagram in the parameter 
space defined by U and a. Depending on the param- 
eters value, different kinds of bound states are formed. 
We studied both their properties and the transition from 
one bipolaron type to another. 



At small electron-phonon coupling two electrons form 
a extremely weakly bound state when U = 0. In this 
regime the binding energy decreases fast with increas- 
ing the Coulomb repulsion and increases with increasing 
phonon frequency. 

For larger couplings, depending on the value of U, the 
phonon-induced attraction may result in the formation 
of a strongly bound SO bipolaron or of a weakly bound 
SI bipolaron. The SO bipolaron forms at small values 
of U, for couplings larger than a c ~ 1. It is an on-site 
state and its properties are well described by the strong 
coupling perturbation theory. With increasing U, around 
the value given by the strong coupling perturbation the- 
ory (U — 2E P ), the SO bipolaron transforms sharply 
(but continuously) into a weakly bound state with an 
enhanced probability to have the electrons on nearest- 
neighbor sites, called SI bipolaron. The SI bipolaron is 
a large state spread over many lattice sites for a in the 
intermediate regime (which corresponds to polaron tran- 
sition region) , and becomes nearest-neighbor localized at 
large a. The unrenormalized exchange energy which wins 
over the reduced polaron effective hopping is responsible 
for binding the SI bipolaron. The binding energy of the 
SI bipolaron and the critical Coulomb repulsion U where 
the bipolaron state disappears are smaller than the val- 
ues obtained in the strong coupling perturbation theory 
(U = 4E P ). 

We found that the ground state always has s-wave sym- 
metry. In the intermediate and strong electron-phonon 
coupling regime, for values of U which correspond to the 
SI ground state, an excited d-wave stable state also ex- 
ists. This state is spatially large for intermediate a and 
nearest- neighbor localized for large a, similar to the cor- 
responding s-wave ground state. The excitation energy 
(difference between the (i-wave and s-wave SI states) is 
larger at intermediate coupling and goes exponentially to 
zero when a is increased. 
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